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We propose a method to analytically show the possibility for the appearance of a maximum in 
the signal-to-noise ratio in nonpotential systems. We apply our results to the FitzHugh-Nagumo 
model under a periodic external forcing, showing that the model exhibits stochastic resonance. The 
procedure that we follow is based on the reduction to a one-dimensional dynamics in the adiabatic 

qq ■ limit, and in the topology of the phase space of the systems under study. Its application to other 

^\ | nonpotential systems is also discussed. 
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I. INTRODUCTION 



Stochastic resonance ( SR ) Jl[- || is a phenomenon in which an enhancement of the response of a non-linear 
system is observed when this system is yielded to an external forcing at some optimized nonzero noise level. Since the 
original proposition of SR as a possible explanation for periodic recurrences in the global climate dynamics, SR has 
become the object of copious theoretical and experimental research in a wide variety of systems in physics, biology 
and chemistry. In all these works the possibility of noise having beneficial effects in the dynamics of nonlinear systems 
H ■ nas been pointed out. The original formulation of the problem, in terms of a bistable system and a periodic forcing 
has been extended to systems under the action of aperiodic forcing W and nondynamical systems W , JrJ . 

In the present work we focus our attention on nonpotential systems. Non-potential systems correspond to systems 
far from equilibrium for which the principle of detailed balance does not hold. There are abundant examples of 
such systems arising from biological, chemical and physical problems. Our contribution in this paper is to develop 
a formalism which allows us to analytically treat a wide class of nonpotential systems among which one can include 
excitable and threshold systems as well as, for example, symmetric double-well models [pd| . In particular, we apply 
. our approach to continue the work undertaken by some authors in studying the stochastic properties of the FHN 
model. This is a well known model with wide application in the field of neuronal research |l3| , || . Apart from several 
numerical simulations done on this subject, Collins et al. || have carried out some analytical work on this matter in 
O ■ the area of aperiodic stochastic resonance . Some experimental research has been performed to show the existence of 
SR in this model |^4j . The results obtained were compared to the predictions of the theory of SR in nondynamical 
systems |^0| . Our scheme allows to analytically approach this problem in a simple way by using a generalization of 
the kinetic equations approach used in the case of potential systems ( see Q , [Q , ) . 

All of the aforementioned models have a common feature; their dynamics exhibit three fixed points: an unstable 
point between two stable points. This feature established some resemblance between the processes described by these 
models and the hopping through a potential barrier. There are a great variety of systems that contains these features. 
The FitzHugh-Nagumo model, in its bistable regime, belongs to them. It is worth pointing out that this is not the 
regime in which this model is used to model neural activity. In this context the FHN model is taken to be in the 
excitable regime where only one globally attractor exists. As it is pointed out by Wiesenfeld et al. (6) , a simple model 
of excitable system consists, among other things, of a threshold or potential barrier. Our theory provides a way to 
compute the escape rates from the attractors of a type of two-dimensional non-potential systems, and therefore it 
furnishes us with the main ingredient to apply the theory developed by Wiesenfeld et al.. In this fact you can find the 
relevance of our work to the field of excitable systems. Another example which fit the characteristics we are asking 
for is the class of symmetric double well systems [|ll| . The Sel'kov model jl6| for autocatalytic systems gathers these 
features, too. 

This paper is organized as follows. In Section II we precisely define the range of applicability of our theory. We 
study the dynamics of the fluctuations and compute the kinetic equations. In Section III we introduce the FitzHugh- 
Nagumo model. We analyze its stochastic properties and show the existence of stochastic resonance. Finally in 
Section IV we discuss our main results. 
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II. DYNAMICS OF THE FLUCTUATIONS 



In this work we study the stochastic properties of a class of two-dimensional nonpotential noisy dynamical systems. 
These system may be characterized by the peculiar topology of the phase space of their corresponding deterministic 
underlying versions. For the case we are concerned the dynamics is characterized by the presence of three aligned 
fixed points; an unstable point between two stable points. An example of this kind was given by Maier and Stein |Tl| ] 
in the context of the escape problem. They studied a system with a symmetric phase space consisting of an hyperbolic 
point between two sinks whose attraction basins were separate by the separatrix of the hyperbolic point. 

To begin with consider a general two-dimensional noisy dynamical system |L2| 

^ = u{x) +$■&), (1) 

where x is the vector whose components are state variables, the field u(x) is the drift, g is the noise matrix and 
is a gaussian white noise of zero mean and correlation function given by 

(&(W)> = 2D5ij5(t-l/), (2) 
with D being the noise intensity. For the case discussed in the components of the drift are 

ui(x) = P 3 {x)+ax m y + b, (3) 



u 2 (x) = cx-dy + e, (4) 

where P$(x) is a third order polynomial, and m is an integer such that < m < 2. From the second of these equations, 
it is easy to check that all the fixed points are aligned. By equating eq. (|^) to zero, one obtains a third order equation 
with three real roots for the proper values of the parameters. 

The corresponding Fokker-Planck equation for the probability density, p(x, t), is 

^ = V-(-u(x)p + V-(Bp)), (5) 

where D = Dg ■ g is the diffusion tensor. 

Let us now assume that the system is potential. In such a case it is possible to write the Fokker-Planck equation 
as a continuity equation: 

g-V./, (0, 

where J is the diffusion current given by 

J = -De- l, / D V -e^ D . (7) 
To obtain this expression, we have defined U and p ( a generalized chemical potential) as follows 

U = - [ u-dx, (8) 



e^ D = <jpe u l D . (9) 

For a potential system, the function U as defined in (j^) is simply the potential energy. In the nonpotential case, 
however, the value of U will depend on the path of integration we choose and, in general, we can not achieve eqs. (|^) 
and®. 

Now consider the weak noise limit and think about some general characteristics of the probability distribution. If 
we let the system evolve during a sufficiently long time the probability distribution will have two maxima at the two 
stable fixed points ( SFP ) of the deterministic dynamics and a minimum in the unstable fixed point ( USFP ) jl7) . 
On the other hand, in the weak noise limit the probability distribution will be very narrow around the line on which 
the fixed point lies. Thus we can assume that the fluctuations run over this line and, therefore, their dynamics are 
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practically one-dimensional. This approximation can be justified accounting for the assumption of low noise level and 
the adiabatic limit. As well as in the one-dimensional double- well problem the adiabatic hypothesis implies that the 
representative point of the system is, in this long time limit, in one of the two wells 0], 0], in our two dimensional 
problem implies that the dynamics will be restricted to be on the nullcline. On the other hand, as it was found by 
Maier and Stein in [Tlj] , the distance to the nullcline, in the limit of weak noise intensity, is normally distributed 
with variance equal to the square root of noise level. Therefore, the fluctuations will run in a very narrow region 
around the nullcline. Consequently, although the whole system is nonpotential, we can reduce the dynamics being 
one-dimensional whose potential is given by the function U taking as the integration path the null cline. For the case 
discussed previously from eqs. (|^), (Q) and (||), under these circumstances ,one then obtains the drift 

u x {x) = P 3 (x) + ^x m+1 +b+^p. (10) 
corresponding to the dynamics of the fluctuations in the weak noise limit and 

u = - f dt( p 3 (t) + ^t m+1 + b + ^H, (ii) 

Our next step will be to discretize the dynamics on the null cline . In particular we will obtain the kinetic equations. 
To this end we define the populations n + ( n_ ) as the population on the right ( left ) of the USFP |is|| . 

p(x,t)dx (12) 

S[+] 



p(x,t)dx (13) 

s[-] 

where S[+] ( S[— ] ) is the portion of the phase space on the right ( left ) of a line orthogonal to the line which contains 
the fixed points and passing through the USFP . 

In the adiabatic limit we can assume that the population is strongly concentrated in a small region around the 
SFP, as suggested by the picture of the probability density that we have profiled when the maxima in this long time 
limit is very high. This corresponds, in the assumption of one-dimensional fluctuations dynamics, having a bistable 
potential with two deep minima at the SFP and maximum at the USFP, or cquivalcntly a high potential barrier. 

In the present context we understand by adiabatic aproximation a long time limit in which all the system has 
arrived to a quasistationary state such that the probability of the system to be in a state different from the stationary 
stable states is practically zero. 

So in this limit we can suppose that the system reaches a quasi-stationary state in which a quasi-stationary diffusion 
current is established. In addition, is assumed to be uniform between the two maxima of the probability density and, 
in the weak noise limit this current is concentrated in the line joining the three fixed points without loss of generality, 
the system can be taken to lie in the axis y = 0. 

J(x,t) = J(x,t)6(y) = J(t)S(y)(9(x - x+) - 9(x - (14) 

where x+ ( x_ ) is the coordinate of the fixed point on the right ( left ) of the USFP. 
The kinetic equation for n + is given by 

dn+ f df> -dx = - / V • Jdx. (15) 



dt 7s[+] dt j S [ + ] 
By using the divergence theorem and the assumptions about the form of the diffusion current, we have 

dfi /*+°° 

—± = J Ji(x,t).6(y)dy = A(x,t), (16) 
and, proceeding in the same way for n_, we obtain 

— = -Ji(x,t) (17) 
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In addition, due to the height of the maxima in the probability density and the weakness of the noise, we can also 
consider that equilibrium in each side of the USFP is reached independently, thus, the generalized chemical potential 
is given by 

H(x,t) = {p(x + ,t)d(x - x) + fj,(x^,t)9(x - x )}S(y), (18) 

where xo is the coordinate of the USFP. The tensorial character of the generalized potential has been removed due 
to the dynamics reduction to one dimension. In equation ( |l8| ) fi = fin has been defined. By using equation ([lj) in 
(0) we obtain 

V(x, t) = V + e-( u - u +V D 9(x -x) + V-e-( u - u -V D 9(x - x )}5(y), (19) 

where U corresponds to the integral over the adequate path of (|l^) , U+ and U- are its values at the SFP ( its minima 
), * = gup and #± = g n p(x±, t). 

In order to obtain the expression for the quasi-stationary current J\{i) we follow the same procedure as in [ EJ . 
From the definition of the probability current and the adiabatic hypothesis we have 

,h(t)(9(x - x+) - 6(x - x_)) = -De- u ' D d x e^ D , (20) 

where a diagonal diffusion tensor is assumed. Integrating now over x and taking into account that, due to the height 
of the barrier, the main contribution to these integrals is around the maximum of the potential xq, one obtains 

J l{ t) = _ jD (M)l/2 e -C/o/C (e M + /S_ e M-/C )j (21) 

2nD 

where U = 4-^- \ xq , Uq = U(xq) and [tT are the values of the generalized chemical potential at SFP. 



u dK^J^JJ' " \ >C_ ' r~ 

By using eq. ( |19| ) we can rewrite (|2lj) in the following way 



, hi t) = D (Bdy/2^_ e ~(u -u-)/D _ 9 -(uo-u+yD^ (22) 



where U+ and C/_ are the values of the potential at the SFP. On the other hand, by using eq. (19) in eq. ( |12[ ) one 
obtains the following relations 

*- = <-S>" j "- < 24 » 

These expressions allow us to write the kinetic equations for the populations n + and rt_ in the form 

dn + dri- 
~~dT ~ dt 

where the kinetic coefficients are given by 



= - K + n + , (25) 



«t = ^(I^X) 1 ^"^-^, (26) 
With this result we have obtained of the kinetic equations for the nonpotential system. 

III. THE FITZHUGH-NAGUMO MODEL 

In this section, we will apply the results of the previous section to the study of the FitzHugh-Nagumo ( FHN 
) neural model @, UJ- |lj. This model is a variant of the Hodgkin-Huxley model |pl| - |p3| which accounts for 
the essentials of the regenerative firing mechanisms in nerve cells. The FHN equations correspond to an excitable 
threshold model but, as will be seen briefly, due to their cubic non-linearity, they exhibit the characteristic behavior 
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of a bistable system. Our main objective is to show analytically the appearance of SR in this model under a periodic 
external forcing. 

The non-dimensional equations of the FHN model are [pi] 

± =v(a-v)(v-l)-w + I a , (27) 
at 



— = bv - 7U>. (28) 
at 

where < a < 1 is essentially the threshold value, b and 7 are positive constants and I a is the applied current. For 
the sake of simplicity, and without loss of generality, we will take I a = 0. The drift field for this model is given by 

ui(v,w) = v(a — v)(v — 1) — w (29) 



U2(v,w) = bv — jw (30) 

As can be seen from eq. ( pj| ) the null cline of the deterministic dynamics of this equations is the line v = ^w. By 
substitution on the right hand side of equation ((2?]) we find the following equation for steady states 

via - v)(v - 1) - -v = 0. (31) 

7 

This is a third order equation, which for certain values of the parameters has three roots ( see Figure 1 ). Among 
these three fixed points two are stable: F_ and F + , and one unstable: F , situated between the other two (^l). Thus, 
in this case, this system fulfills the conditions in order to our theory to be applied. 

The function U defined in Section II has to be integrated in this case over line v = ^u>. Performing this integration 
one obtains 

U = — - <£±^ + <Z±Mj. (32) 
4 3 2 

On the other hand, when this system is in a noisy environment, in the limit of weak noise, we can approximate the 
dynamics of the fluctuations by the one-dimensional Langevin equation 

^ = v (a-v)(v-l)-^v + at), (33) 

that is, the fluctuations run along the line v = Zw. As can be easily checked, U is the potential for the deterministic 
part of this equation. The two SFP of this one-dimensional dynamics are the two minima of ( |32"|) and the USFP is 
its maximum. Collins et al. || have arrived previously to similar conclusions by another approach in the context of 
the study of aperiodic stochastic resonance. 

Fig. 2 shows the asymmetric shape of U. Before going any further, it would be interesting to summarize what 
this picture has to say to us about the physics of the problem ||. In the FHN model there is a fast variable, 
v(t), and a recovery-like variable, w(t). After the barrier threshold has been crossed, i.e. the system has gone to an 
"excited" state, the system returns ( even in the deterministic case ) to the "rest" state. As can be seen in Fig. 2, 
there is one of the stable states for which the potential is larger than for the other stable state. Therefore, there is a 
more stable state, which corresponds to the rest state to which the system after some time , under the action of the 
noise, returns. In our scheme the presence of noise is necessary in order to return to the rest state, because of the 
elimenation of w(t). 

In order to show how this scheme can account for the existence of stochastic resonance in the FHN model, we 
will suppose that the system is under the action of periodic forcing. For simplicity's sake, we will assume that the 
parameter a is a periodic function: a = ao(l + eosinw s t) where £q is a small parameter and ao(l + a) < 1. 

To take a as an oscillatory factor implies that the positions of Fq ( the USFP ) and F+ ( one of the SFP ) as well 
as the values of the potential at these points become oscillatory functions, too. The position of F_ (the other SFP 
) remains constant. Let Vq, V— and w + be the v-coordinate of the maximum and the minima, respectively, of the 
potential; one has 

U+ =U(v+) = C++ v+<t). (34) 
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To compute the moments and the power spectrum, we assume that, in the limit of weak noise, the probability 
density ( in one dimension ) can be written as (2) 

p(v,t) = n + (t)6 v ^ v+ + ra_ (*)£„,„_, (35) 

where n + and n_ come from the kinetic equations (5a). The formal solution to these equations is found to be 



n±(t) = g- 1 (t)(n±(t Q )g(t Q )+ f K T (t')g(t')dt'), (36) 

J t 

with 

g(t) = J (K + + K_)dt'. (37) 

In order to calculate n± we perform a Taylor expansion of the transition rates up to the first order in respect to 
the parameter e(t). 

K± = + aftfiosinuot (38) 
where <fio has been defined as (f>o = sq/D and and af are given by 

p-io/D 

a± = (39) 
' " !/fl -4e^/ fl (r,o-r, ± ), (40) 



2vr 

with being the zero order contribution of \Uq\U±, £0 and f?o the zero and first order contribution of U at the USFP 
and £j- and rj± the zero and first order contribution of U at the SFPs. Its first order contribution can be neglected in 
the limit of weak noise. By introducing these Taylor expansions in equations ([36|) and (^), we get the populations 
up to the first order in the parameter 4>q. 



n±(t) - e- 00 ^ 10 ^^^),^ - -SL - fa-U—Lcosfat) + foal 



Y'U ^u^u„y 1 rU uj „ 2M/9 



ai cos(u; t -HE) «i , cos(tJot + $) 

~ ^,2x,,2U/2 ) + — ^ + ^0— cos(wot)) - foa+ + 

(41) 



u (og + wg)Va ' « V ru Wo v u " ru 1 (ag + wg)Va 
ai sin(uJot + $) 



«;„ (ag + a, 2 ) 1 /*' 



where $ = arctg(ao/u>o) and c*o = Q<o + Q!^". The quantity n±(t\vt , to) is the conditional probability that v(t) is 
in the + state at time t given that the state at time to was i>j D . From this equation it is possible to compute the 
statistical properties of the process v(t). Of particular interest to our purposes is to find its autocorrelation function 
which is given by (Wiesenfeld) 

(v(t)v(t + t)) = I (v(t)v(t + r)\vt ,to). (42) 

The conditional correlation function is given by 

(v(t)v(t + r)\v ta ,t Q ) = v+(t + T)v + (t)n + (t + T\v + ,t)n + (t\v to ,to) + 
v+(t + r)v-(t)n + (t + T\v_,t)n_(t\v to ,t ) + 
V-(t + T)v + (t)n-(t + T\v + ,t)n + (t\v to ,to) + 

v_ (t + t)v- (t)n_ (t + t\v+, (t\v to , to). (43) 

Let us make some considerations which allows us to simplify the computation of the autocorrelation function. It 
is clear from equation (43) that the Fourier transform of the autocorrelation function will depend on t as well as 
in the frequency. This dependency is avoided by taking its average over the period of the external forcing [Q. The 
autocorrelation function is to be computed up to the second order in parameter O ~ -D _1 .Thus, in the limit of the 
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weak noise D" 1 , can be neglected in comparision to D~ 2 . Therefore, the only contribution of the first order term of 
v± to the autocorrelation function comes from its product with the zero order term of the product of n's. But, on 
doing the average this term vanishes. So, finally, and taking into account that u_ = 0, we arrive to 

(v(t)v(t + r)\v to ,t ) = \\n + (t + T\v + ,t)n + (t\v to ,t ), (44) 

where the overlined indicates an average over t and A+ is the position of F + up to order zero. From equations (43) 
and ( [44| ) taking the average and the limit to —> —00, we finally obtain the following expression for the autocorrelation 
function. 



(v(t)v(t + r)) = A^) 2 + e^l(^)(l-^) + 

A +^o{o e a 0\ — > - )(-2-j 2 + 

Oo~a{ax ai 2 cos(u>ot) 

\ a o{ — j j— —. — 2 + 

a luq Uo + w o 

(— — + 1 )sin(uj T)) + 
a u>o a 

cos{uj t) f a Q a 2 a^a^ai | 1 ( _ 2 1 «i 2 

2 1 2 V 2 + o\ 1 + o > IS' 

With this results we can now compute the averaged power spectrum given by 



(45) 



5(0) = ^ / S(n,t)dt = I {v(t)v(t + T )}e~ tUT dT, (46) 

where the last equality follows from the commutative character of averaging and Fourier transform. We obtain, after 
Fourier transform (50) 

SW) = ^6(Q) + 2-tA s Oo-(1 - ^) + 2 7 r^±4(i^ - 
v ' + a Q v ' al + il 2 v a ; a 2 , + O 2 v 

(Q °' )2 f ai + 1(^5! + k^mn - co). (47) 

Uq 2 a 2 cj 

In this expression the fraction of the total power in the broadband noisy part of the spectrum, which usually is a 
small fraction of the total power, has been neglected 0. Note that in the power spectrum there is a term proportional 
to 5(£l). This is due to the asymmetry of the potential that originates a mean probability current between one 
stable state and the more stable one is more stable than the other one. This term has not been obtained in previous 
approaches to this problem |{J , flo|| 

From equation (47) the signal to noise ratio, R, can be obtained as a function of the noise level D, by making 
Q = Wo- 

R = ^ ^ . (48) 

1 _ 2a. 

We have plotted R in Figure 3 as a function of the noise level. It exhibits a maximum for a non-zero noise level 
and therefore this model shows stochastic resonance. 



IV. CONCLUSIONS 



In this paper we have proposed a method to analyze the possibility for the appearance of SR in nonpotential 
systems. In particular we have treated the FHN model. A lot of numerical work on this subject has been performed 
but only a few papers have treated this problem from an analytical perspective due to the difficulties inherent to the 
nature of nonpotential system. The work by Collins et al fj| where aperiodic SR is discussed, and the paper by Pei 
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et al. |24|], in which the theoretical framework developed by Gingl et al. ]10| was used to interpret some experimental 
results, are among the ones discussing analytical aspects. 

Far from being specific of this problem, the method we have proposed can be extended to a wide variety of different 
nonpotential systems ranging from threshold systems to some autocatalytic models and symmetric double well systems 
from fields so diverse as chemical kinetics and population dynamics. Although, for the sake of brevity, we do not quote 
our results here we have applied our approach to that of Collins et al. || with a periodic input and we obtain similar 
results in both cases, not the same because in the case of these authors the input is additive and the threshold is 
maintained constant. We have applied our theory to the standard double- well model and we have obtain that this 
model exhibits SR. The result in this case is equivalent to the one obtained in Q for a symmetric quartic potential. 

The scheme we present in this paper allows us to treat potential and the class of aforementioned nonpotential 
systems in an unified way. Moreover, it reproduces the essentials of the physics of the problem as can be seen from 
the obtaining of the refractory current in eq. (52) ]l3| ]. 

Finally, it is worth pointing out that our method, which constitutes essentially an extension of the Kramers rate 
theory to this kind of nonpotential systems, enables one to compute the kinetic coefficients in a simple and direct 
way. 
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FIGURE CAPTIONS 

• Figure 1.- Fixed points of the FitzHugh-Nagumo model for a = 0.5 and ^ = 0.01. Non- 
dimensional variables. 

• Figure 2.- Asymmetric one-dimensional potential as a function of v. Non-dimensional variables. 

• Figure 3.- Signal-to-Noise Ratio as a function of the noise level, D. Non-dimensional variables. 
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